rm(list = ls())
print(getwd())
## Load libraries

library(ggplot2, warn.conflicts = FALSE)
library(data.table, warn.conflicts = FALSE)
library(plm, warn.conflicts = FALSE)
library(gridExtra, warn.conflicts = FALSE)
suppressPackageStartupMessages(library(texreg, warn.conflicts = FALSE))
suppressPackageStartupMessages(library(panelView, warn.conflicts = FALSE))
suppressPackageStartupMessages(library(tidyverse, warn.conflicts = FALSE))

## Load data

load("../generated_data/bes_w1_w21_pacer_w1_w2_merged.Rdata")
load(file = "../generated_data/fe_models_new.Rdata")

setkey(bes_long, id, wave)

## #################
## Outputs for paper

## - Baseline spec = collapse UC and furlough

get_marginal_effects <- function(x = NULL){
  
  fx_mean <- coef(x)
  fx_se <- sqrt(diag(vcov(x)))
  fx_hi <- fx_mean + 1.96 * fx_se
  fx_lo <- fx_mean - 1.96 * fx_se
  
  est_out <- data.frame(term = names(coef(x)),
                        estimate = fx_mean,
                        se = fx_se,
                        conf.high = fx_hi,
                        conf.low = fx_lo,
                        row.names = NULL)
  
  est_out <- tibble(est_out)
  
  return(est_out)
  
}


get_marginal_effects_interaction <- function(var, interaction_levels = levels(out$vote19), x = NULL){
  
  var_coefs_idx <- which(grepl(var,names(x$coefficients)))
  var_coefs_baseline_idx <- var_coefs_idx[1]
  var_ceofs_int_inx <- var_coefs_idx[-1]
  
  fx_mean <- c(x$coefficients[var_coefs_baseline_idx],
               x$coefficients[var_coefs_baseline_idx] + x$coefficients[var_ceofs_int_inx])
  fx_se <- c(sqrt(vcov(x)[var_coefs_baseline_idx, var_coefs_baseline_idx]), 
             sapply(1:length(var_ceofs_int_inx), function(i) sqrt(vcov(x)[var_coefs_baseline_idx, var_coefs_baseline_idx] + vcov(x)[var_ceofs_int_inx[i], var_ceofs_int_inx[i]] + vcov(x)[var_coefs_baseline_idx,var_ceofs_int_inx[i]])))
  
  fx_hi <- fx_mean + 1.96 * fx_se
  fx_lo <- fx_mean - 1.96 * fx_se
  
  est_out <- data.frame(treatment = var,
                        levels = interaction_levels,
                        estimate = fx_mean,
                        se = fx_se,
                        hi = fx_hi,
                        lo = fx_lo,
                        row.names = NULL)
  
  return(est_out)
  
  
  
}

baseline_out_income_guarantee_coefs <- baseline_out_income_guarantee %>% 
  map(~ get_marginal_effects(.x)) %>% 
  bind_rows(.id = "var") %>%
  mutate(term = recode(term, lost_employmentTRUE = "Lost Employment",
                       income_guarantee_treatmentTRUE = "Gov. Support"))

# 1. onelawrich; taxSpendSelf; redistSelf (Aii)

outcome_vars <- c("redistSelf", "taxSpendSelf", "ideoOneLawRich", "ideoBHPSJjobForAll", "reasonForUnemployment", "riskPoverty", "riskUnemployment", "govtHandouts")

baseline_out_income_guarantee_coefs %>% 
  filter(var %in% c("ideoBHPSJjobForAll","taxSpendSelf", "redistSelf") & term == "Gov. Support") %>%
  mutate(var = factor(var, levels = outcome_vars),
         var = recode_factor(var, `ideoBHPSJjobForAll`= "jobsForAll")) %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, y = var)) + 
  geom_point() +
  geom_errorbarh(height = 0.001) + 
  geom_vline(xintercept = 0, linetype = 2) + 
  theme_bw() +
  theme(axis.text = element_text(size = 12),
        text = element_text(size = 12)) + 
  ylab("") + 
  xlab("Treatment effect") +
  xlim(c(-1,1)) + 
  facet_wrap(~term)

ggsave("../generated_images/fig2_fe_baseline_ideo_vars.pdf", height = 5*.75, width = 10)
ggsave("../generated_images/fig2_fe_baseline_ideo_vars.tiff", height = 5*.75, width = 10, dpi = 800)


# Table
ns <- rbind(apply(attr(baseline_out_income_guarantee$redistSelf$model, "index"),2,function(x) length(unique(x))),
            apply(attr(baseline_out_income_guarantee$taxSpendSelf$model, "index"),2,function(x) length(unique(x))),
            apply(attr(baseline_out_income_guarantee$ideoBHPSJjobForAll$model, "index"),2,function(x) length(unique(x))))
n_obs <- c(length(baseline_out_income_guarantee$redistSelf$weights),
           length(baseline_out_income_guarantee$taxSpendSelf$weights),
           length(baseline_out_income_guarantee$ideoBHPSJjobForAll$weights))

sink("../generated_tex/fe_baseline_ideo_vars.tex")
m1 <- baseline_out_income_guarantee$redistSelf 
m2 <- baseline_out_income_guarantee$taxSpendSelf
m3 <- baseline_out_income_guarantee$ideoBHPSJjobForAll
stargazer::stargazer(m1, m2, m3,
                     dep.var.labels.include = T,
                     dep.var.caption = "",
                     multicolumn = FALSE,
                     dep.var.labels = c("$redistSelf$", "$taxSpendSelf$", "$jobsForAll$"),
                     no.space = TRUE,
                     omit.stat = c("all"),
                     covariate.labels = c("$CovidLostEmployment$", "$CovidGovSupport$"),
                     float = FALSE,
                     add.lines = list(c("N respondents", ns[,1]),
                                      c("N waves", ns[,2]),
                                      c("N obs.", n_obs),
                                      c("Ind. FEs", "Yes", "Yes", "Yes"),
                                      c("Wave FEs", "Yes", "Yes", "Yes")))
sink()



# 2. RiskUnemployment and RiskPoverty (Aiii)

baseline_out_income_guarantee_coefs %>% 
  filter(var %in% c("riskPoverty", "riskUnemployment")) %>%
  mutate(var = factor(var, levels = outcome_vars)) %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, y = var)) + 
  geom_point() +
  geom_errorbarh(height = 0.001) + 
  geom_vline(xintercept = 0, linetype = 2) + 
  theme_bw() +
  ylab("") + 
  xlab("Treatment effect") +
  xlim(c(-1,1)) + 
  facet_wrap(~term)

ggsave("../generated_images/fig5_fe_baseline_risk_vars.pdf", height = 2, width = 6)
ggsave("../generated_images/fig5_fe_baseline_risk_vars.tiff", height = 2, width = 6, dpi = 800)

# Table

ns <- rbind(apply(attr(baseline_out_income_guarantee$riskPoverty$model, "index"),2,function(x) length(unique(x))),
            apply(attr(baseline_out_income_guarantee$riskUnemployment$model, "index"),2,function(x) length(unique(x))))
n_obs <- c(length(baseline_out_income_guarantee$riskPoverty$weights),
           length(baseline_out_income_guarantee$riskUnemployment$weights))

sink("../generated_tex/fe_baseline_risk_vars.tex")
m1 <- baseline_out_income_guarantee$riskPoverty 
m2 <- baseline_out_income_guarantee$riskUnemployment
stargazer::stargazer(m1, m2, 
                     dep.var.labels.include = T,
                     dep.var.caption = "",
                     multicolumn = FALSE,
                     dep.var.labels = c("$riskPoverty$", "$riskUnemployment$"),
                     no.space = TRUE,
                     omit.stat = c("all"),
                     covariate.labels = c("$CovidLostEmployment$", "$CovidGovSupport$"),
                     float = FALSE,
                     add.lines = list(c("N respondents", ns[,1]),
                                      c("N waves", ns[,2]),
                                      c("N obs.", n_obs),
                                      c("Ind. FEs", "Yes", "Yes", "Yes"),
                                      c("Wave FEs", "Yes", "Yes", "Yes")))
sink()


# 3. govtHandouts and reasonForUnemployment (Deservingness)

baseline_out_income_guarantee_coefs %>% 
  filter(var %in% c("govtHandouts", "reasonForUnemployment")) %>%
  mutate(var = factor(var, levels = outcome_vars)) %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, y = var)) + 
  geom_point() +
  geom_errorbarh(height = 0.001) + 
  geom_vline(xintercept = 0, linetype = 2) + 
  theme_bw() +
  ylab("") + 
  xlab("Treatment effect") +
  xlim(c(-1,1)) + 
  facet_wrap(~term)

ggsave("../generated_images/fe_baseline_deserving_vars.pdf", height = 4/2, width = 12)



### - Baseline spec = separating UC, furlough, self-employment

baseline_out_income_guarantee_separate_coefs <- baseline_out_income_guarantee_separate %>% 
  map(~ get_marginal_effects(.x)) %>% 
  #map(~ broom::tidy(.x,conf.int = TRUE)) %>% 
  bind_rows(.id = "var") %>%
  mutate(term = recode(term, lost_employmentTRUE = "Lost Employment",
                       income_guarantee_covid_job_treatmentTRUE = "Furlough",
                       income_guarantee_covid_self_treatmentTRUE = "Self-employment",
                       income_guarantee_other_treatmentTRUE = "Universal Credit"))

# 4. onelawrich; taxSpendSelf; redistSelf (Bi)

baseline_out_income_guarantee_separate_coefs %>% 
  filter(var %in% c("ideoBHPSJjobForAll","taxSpendSelf", "redistSelf")) %>%
  mutate(var = factor(var, levels = outcome_vars),
         var = recode_factor(var, `ideoBHPSJjobForAll`= "jobsForAll")) %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, y = var)) + 
  geom_point() +
  geom_errorbarh(height = 0.001) + 
  geom_vline(xintercept = 0, linetype = 2) + 
  theme_bw() +
  ylab("") + 
  xlab("Treatment effect") +
  xlim(c(-1,1)) + 
  facet_wrap(~term, ncol = 2)

ggsave("../generated_images/fe_multi_treatment_ideo_vars.pdf", width = 6, height = 3)

# Table
ns <- rbind(apply(attr(baseline_out_income_guarantee_separate$redistSelf$model, "index"),2,function(x) length(unique(x))),
            apply(attr(baseline_out_income_guarantee_separate$taxSpendSelf$model, "index"),2,function(x) length(unique(x))),
            apply(attr(baseline_out_income_guarantee_separate$ideoBHPSJjobForAll$model, "index"),2,function(x) length(unique(x))))
n_obs <- c(length(baseline_out_income_guarantee_separate$redistSelf$weights),
           length(baseline_out_income_guarantee_separate$taxSpendSelf$weights),
           length(baseline_out_income_guarantee_separate$ideoBHPSJjobForAll$weights))

sink("../generated_tex/fe_multi_treatment_ideo_vars.tex")
m1 <- baseline_out_income_guarantee_separate$redistSelf
m2 <- baseline_out_income_guarantee_separate$taxSpendSelf
m3 <- baseline_out_income_guarantee_separate$ideoBHPSJjobForAll
stargazer::stargazer(m1,m2,m3,
                     dep.var.labels.include = T,
                     dep.var.caption = "",
                     multicolumn = FALSE,
                     dep.var.labels = c("$redistSelf$", "$taxSpendSelf$", "$jobsForAll$"),
                     no.space = TRUE,
                     omit.stat = c("all"),
                     covariate.labels = c("$CovidLostEmployment$", "$CovidFurlough$", "$CovidSelfEmployment$", "$CovidUniversalCredit$"),
                     float = FALSE,
                     add.lines = list(c("N respondents", ns[,1]),
                                      c("N waves", ns[,2]),
                                      c("N obs.", n_obs),
                                      c("Ind. FEs", "Yes", "Yes", "Yes"),
                                      c("Wave FEs", "Yes", "Yes", "Yes")))
sink()


# 5. govtHandouts and reasonForUnemployment (Deservingness)

baseline_out_income_guarantee_separate_coefs %>% 
  filter(var %in% c("govtHandouts", "reasonForUnemployment") ) %>%
  mutate(var = factor(var, levels = outcome_vars)) %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, y = var)) + 
  geom_point() +
  geom_errorbarh(height = 0.001) + 
  geom_vline(xintercept = 0, linetype = 2) + 
  theme_bw() +
  ylab("") + 
  xlab("Treatment effect") +
  xlim(c(-1,1)) + 
  facet_wrap(~term, ncol = 2)

ggsave("../generated_images/fig7_fe_multi_treatment_deserving_vars.pdf", width = 6, height = 3)
ggsave("../generated_images/fig7_fe_multi_treatment_deserving_vars.tiff", width = 6, height = 3, dpi = 800)


# Table
ns <- rbind(apply(attr(baseline_out_income_guarantee_separate$govtHandouts$model, "index"),2,function(x) length(unique(x))),
            apply(attr(baseline_out_income_guarantee_separate$reasonForUnemployment$model, "index"),2,function(x) length(unique(x))))
n_obs <- c(length(baseline_out_income_guarantee_separate$govtHandouts$weights),
           length(baseline_out_income_guarantee_separate$reasonForUnemployment$weights))

sink("../generated_tex/fe_multi_treatment_deserving_vars.tex")
m1 <- baseline_out_income_guarantee_separate$govtHandouts
m2 <- baseline_out_income_guarantee_separate$reasonForUnemployment
stargazer::stargazer(m1, m2, 
                     dep.var.labels.include = T,
                     dep.var.caption = "",
                     multicolumn = FALSE,
                     dep.var.labels = c("$govtHandouts$", "$reasonForUnemployment$"),
                     no.space = TRUE,
                     omit.stat = c("all"),
                     covariate.labels = c("$CovidLostEmployment$", "$CovidFurlough$", "$CovidSelfEmployment$", "$CovidUniversalCredit$"),
                     float = FALSE,
                     add.lines = list(c("N respondents", ns[,1]),
                                      c("N waves", ns[,2]),
                                      c("N obs.", n_obs),
                                      c("Ind. FEs", "Yes", "Yes", "Yes"),
                                      c("Wave FEs", "Yes", "Yes", "Yes")))
sink()


## - party interaction spec = separating UC, furlough, self-employment (might be slicing the data too thin)

party_out_income_guarantee_separate_coefs <- party_out_income_guarantee_separate %>%
  map(~ lapply(c("income_guarantee_covid_job_treatment", "income_guarantee_covid_self_treatment", "income_guarantee_other_treatment"), 
               function(v) 
                 get_marginal_effects_interaction(v, interaction_levels = levels(out$vote19), x = .x)) %>%
        bind_rows) %>% 
  bind_rows(.id = "var") %>%
  tibble() %>%
  mutate(treatment = recode(treatment, "income_guarantee_covid_job_treatment" = "Furlough",
                            "income_guarantee_covid_self_treatment" = "Self-employment",
                            "income_guarantee_other_treatment" = "Universal Credit"))


# 6. Onelawrich; TaxSpendSelf; RedistSelf; ReasonforUnemployment (Bi)

party_out_income_guarantee_separate_coefs %>%
  filter(var %in% c("ideoBHPSJjobForAll","taxSpendSelf", "redistSelf")) %>%
  mutate(var = factor(var, levels = outcome_vars),
         var = recode_factor(var, `ideoBHPSJjobForAll`= "jobsForAll")) %>%
  ggplot(aes(y = var, x = estimate, xmin = lo, xmax = hi, col = levels, group = levels, linetype = levels)) +
  geom_point(position = position_dodge(width = .5)) +
  geom_errorbarh(height = 0.0001, position = position_dodge(width = .5)) + 
  facet_wrap(~ treatment, ncol = 4) + 
  geom_vline(xintercept = 0, linetype = 1) + 
  xlim(c(-1,1)) + 
  theme_bw() + 
  scale_color_manual("Party", values = c("black", "darkgrey", "lightgrey")) + 
  xlab("Treatment effect") + 
  labs(color  = "Party", linetype = "Party") +
  ylab("")

ggsave("../generated_images/fig3_fe_by_party_by_treat.pdf", width = 6, height = 3)
ggsave("../generated_images/fig3_fe_by_party_by_treat.tiff", width = 6, height = 3, dpi = 800)


# Table
ns <- rbind(apply(attr(party_out_income_guarantee_separate$redistSelf$model, "index"),2,function(x) length(unique(x))),
            apply(attr(party_out_income_guarantee_separate$taxSpendSelf$model, "index"),2,function(x) length(unique(x))),
            apply(attr(party_out_income_guarantee_separate$ideoBHPSJjobForAll$model, "index"),2,function(x) length(unique(x))))
n_obs <- c(length(party_out_income_guarantee_separate$redistSelf$weights),
           length(party_out_income_guarantee_separate$taxSpendSelf$weights),
           length(party_out_income_guarantee_separate$ideoBHPSJjobForAll$weights))

sink("../generated_tex/fe_by_party_by_treat.tex")
m1 <- party_out_income_guarantee_separate$redistSelf
m2 <- party_out_income_guarantee_separate$taxSpendSelf
m3 <- party_out_income_guarantee_separate$ideoBHPSJjobForAll
stargazer::stargazer(m1,m2,m3,
                     dep.var.labels.include = T,
                     dep.var.caption = "",
                     multicolumn = FALSE,
                     dep.var.labels = c("$redistSelf$", "$taxSpendSelf$", "$jobsForAll$"),
                     no.space = TRUE,
                     omit.stat = c("all"),
                     omit = c("^vote19Labour:as.factor",
                              "^vote19Other:as.factor"),
                     covariate.labels = c("$CovidLostEmployment$", 
                                          "$CovidFurlough$", 
                                          "$CovidSelfEmployment$", 
                                          "$CovidUniversalCredit$",
                                          "$CovidFurlough * Labour$", 
                                          "$CovidFurlough * Other$", 
                                          "$CovidSelfEmployment * Labour$", 
                                          "$CovidSelfEmployment * Other$", 
                                          "$CovidUniversalCredit * Labour$", 
                                          "$CovidUniversalCredit * Other$"),
                     float = FALSE,
                     add.lines = list(c("N respondents", ns[,1]),
                                      c("N waves", ns[,2]),
                                      c("N obs.", n_obs),
                                      c("Ind. FEs", "Yes", "Yes", "Yes"),
                                      c("Wave FEs", "Yes", "Yes", "Yes")))
sink()


# 7. govtHandouts and reasonForUnemployment (Deservingness)

party_out_income_guarantee_separate_coefs %>%
  filter(var %in% c("govtHandouts", "reasonForUnemployment")) %>%
  mutate(var = factor(var, levels = outcome_vars)) %>%
  ggplot(aes(y = var, x = estimate, xmin = lo, xmax = hi, col = levels, group = levels)) +
  geom_point(position = position_dodge(width = .5)) +
  geom_errorbarh(height = 0.0001, position = position_dodge(width = .5)) + 
  facet_wrap(~ treatment, ncol = 4) + 
  geom_vline(xintercept = 0, linetype = 1) + 
  xlim(c(-1,1)) + 
  theme_bw() + 
  scale_color_manual("Party", values = c("blue", "red", "gray")) + 
  xlab("Treatment effect") + 
  ylab("")

ggsave("../generated_images/fe_by_party_by_treat_deserving.pdf", width = 15, height = 5*.5)


## - prior position interaction spec = separating UC, furlough, self-employment

prior_out_income_guarantee_separate %>%
  map(~ lapply(c("income_guarantee_covid_job_treatment", "income_guarantee_covid_self_treatment", "income_guarantee_other_treatment","lost_employment"), 
               function(v) 
                 get_marginal_effects_interaction(v, interaction_levels = c("Left","Mid", "Right"), x = .x)) %>%
        bind_rows) %>% 
  bind_rows(.id = "var") %>%
  tibble() %>%
  mutate(treatment = recode(treatment, "income_guarantee_covid_job_treatment" = "Furlough",
                            "income_guarantee_covid_self_treatment" = "Self-employment",
                            "income_guarantee_other_treatment" = "Universal Credit",
                            "lost_employment" = "Lost Employment")) %>%
  ggplot(aes(x = estimate, xmin = lo, xmax = hi, y = var, col = levels, group = levels)) + geom_pointrange(position = position_dodge(.2)) + facet_wrap(~treatment) + theme_bw() + geom_vline(xintercept = 0, lty = 2)   + 
  ylab("") + 
  xlab("Conditional treatment effect") + 
  scale_color_manual("Prior attitude", values = c("orange", "darkgreen", "brown"))

ggsave("../generated_images/fe_by_prior_by_treat.pdf", width = 6, height = 4)

## - prior position interaction spec, w14 onwards = separating UC, furlough, self-employment

prior_out_income_guarantee_separate14 %>%
  map(~ lapply(c("income_guarantee_covid_job_treatment", "income_guarantee_covid_self_treatment", "income_guarantee_other_treatment","lost_employment"), 
               function(v) 
                 get_marginal_effects_interaction(v, interaction_levels = c("Left","Mid", "Right"), x = .x)) %>%
        bind_rows) %>% 
  bind_rows(.id = "var") %>%
  tibble() %>%
  mutate(treatment = recode(treatment, "income_guarantee_covid_job_treatment" = "Furlough",
                            "income_guarantee_covid_self_treatment" = "Self-employment",
                            "income_guarantee_other_treatment" = "Universal Credit",
                            "lost_employment" = "Lost Employment")) %>%
  ggplot(aes(x = estimate, xmin = lo, xmax = hi, y = var, col = levels, group = levels)) + geom_pointrange(position = position_dodge(.2)) + facet_wrap(~treatment) + theme_bw() + geom_vline(xintercept = 0, lty = 2)   + 
  ylab("") + 
  xlab("Conditional treatment effect") + 
  scale_color_manual("Prior attitude", values = c("orange", "darkgreen", "brown"))

ggsave("../generated_images/fe_by_prior_by_treat_w14onwards.pdf", width = 6, height = 4)

## - prior position interaction spec, most recent prior position = separating UC, furlough, self-employment

prior_out_income_guarantee_separate_recent %>%
  map(~ lapply(c("income_guarantee_covid_job_treatment", "income_guarantee_covid_self_treatment", "income_guarantee_other_treatment","lost_employment"), 
               function(v) 
                 get_marginal_effects_interaction(v, interaction_levels = c("Left","Mid", "Right"), x = .x)) %>%
        bind_rows) %>% 
  bind_rows(.id = "var") %>%
  tibble() %>%
  mutate(treatment = recode(treatment, "income_guarantee_covid_job_treatment" = "Furlough",
                            "income_guarantee_covid_self_treatment" = "Self-employment",
                            "income_guarantee_other_treatment" = "Universal Credit",
                            "lost_employment" = "Lost Employment")) %>%
  ggplot(aes(x = estimate, xmin = lo, xmax = hi, y = var, col = levels, group = levels)) + geom_pointrange(position = position_dodge(.2)) + facet_wrap(~treatment) + theme_bw() + geom_vline(xintercept = 0, lty = 2)   + 
  ylab("") + 
  xlab("Conditional treatment effect") + 
  scale_color_manual("Prior attitude", values = c("orange", "darkgreen", "brown"))

ggsave("../generated_images/fe_by_prior_by_treat_most_recent.pdf", width = 6, height = 4)

## - "don't know" prior position interaction spec = separating UC, furlough, self-employment

dk_out_income_guarantee_separate %>%
  map(~ lapply(c("income_guarantee_covid_job_treatment", "income_guarantee_covid_self_treatment", "income_guarantee_other_treatment","lost_employment"), 
               function(v) 
                 get_marginal_effects_interaction(v, interaction_levels = c("FALSE", "TRUE"), x = .x)) %>%
        bind_rows) %>% 
  bind_rows(.id = "var") %>%
  tibble() %>%
  mutate(treatment = recode(treatment, "income_guarantee_covid_job_treatment" = "Furlough",
                            "income_guarantee_covid_self_treatment" = "Self-employment",
                            "income_guarantee_other_treatment" = "Universal Credit",
                            "lost_employment" = "Lost Employment")) %>%
  ggplot(aes(x = estimate, xmin = lo, xmax = hi, y = var, col = levels, group = levels)) + geom_pointrange(position = position_dodge(.2)) + facet_wrap(~treatment) + theme_bw() + geom_vline(xintercept = 0, lty = 2)   + 
  ylab("") + 
  xlab("Conditional treatment effect") + 
  scale_color_manual("'Don't know' respondent", values = c("orange", "darkgreen"))

ggsave("../generated_images/fe_by_prior_dk_by_treat.pdf", width = 7.5, height = 5)

## - political attention interaction spec = separating UC, furlough, self-employment

attention_out_income_guarantee_separate %>%
  map(~ lapply(c("income_guarantee_covid_job_treatment", "income_guarantee_covid_self_treatment", "income_guarantee_other_treatment","lost_employment"), 
               function(v) 
                 get_marginal_effects_interaction(v, interaction_levels = c("Low", "High"), x = .x)) %>%
        bind_rows) %>% 
  bind_rows(.id = "var") %>%
  tibble() %>%
  mutate(treatment = recode(treatment, "income_guarantee_covid_job_treatment" = "Furlough",
                            "income_guarantee_covid_self_treatment" = "Self-employment",
                            "income_guarantee_other_treatment" = "Universal Credit",
                            "lost_employment" = "Lost Employment")) %>%
  ggplot(aes(x = estimate, xmin = lo, xmax = hi, y = var, col = levels, group = levels)) + geom_pointrange(position = position_dodge(.2)) + facet_wrap(~treatment) + theme_bw() + geom_vline(xintercept = 0, lty = 2)   + 
  ylab("") + 
  xlab("Conditional treatment effect") + 
  scale_color_manual("Political attention", values = c("orange", "darkgreen"))

ggsave("../generated_images/fe_by_attention_by_treat.pdf", width = 7.5, height = 5)

## - prior government support interaction spec = separating UC, furlough, self-employment

prior_gov_support_out_income_guarantee_separate %>%
  map(~ lapply(c("income_guarantee_covid_job_treatment", "income_guarantee_covid_self_treatment", "income_guarantee_other_treatment","lost_employment"), 
               function(v) 
                 get_marginal_effects_interaction(v, interaction_levels = c("No", "Yes"), x = .x)) %>%
        bind_rows) %>% 
  bind_rows(.id = "var") %>%
  tibble() %>%
  mutate(treatment = recode(treatment, "income_guarantee_covid_job_treatment" = "Furlough",
                            "income_guarantee_covid_self_treatment" = "Self-employment",
                            "income_guarantee_other_treatment" = "Universal Credit",
                            "lost_employment" = "Lost Employment")) %>%
  ggplot(aes(x = estimate, xmin = lo, xmax = hi, y = var, col = levels, group = levels)) + geom_pointrange(position = position_dodge(.2)) + facet_wrap(~treatment) + theme_bw() + geom_vline(xintercept = 0, lty = 2)   + 
  ylab("") + 
  xlab("Conditional treatment effect") + 
  scale_color_manual("Prior government support", values = c("orange", "darkgreen"))

#lapply(prior_gov_support_out_income_guarantee_separate, summary)

ggsave("../generated_images/fe_by_prior_gov_support_by_treat.pdf", width = 7.5, height = 5)


## - interacting lost_employment and universal credit receipt

# Table

sink("../generated_tex/fe_multi_treatment_employment_uc_interaction.tex")
m1 <- baseline_out_income_guarantee_separate_lost_employment_interaction$redistSelf
m2 <- baseline_out_income_guarantee_separate_lost_employment_interaction$taxSpendSelf
m3 <- baseline_out_income_guarantee_separate_lost_employment_interaction$ideoBHPSJjobForAll
m4 <- baseline_out_income_guarantee_separate_lost_employment_interaction$riskUnemployment
m5 <- baseline_out_income_guarantee_separate_lost_employment_interaction$riskPoverty

stargazer::stargazer(m1, m2, m3, m4, m5,
                     dep.var.labels.include = T,
                     dep.var.caption = "",
                     multicolumn = FALSE,
                     dep.var.labels = c("$redistSelf$", "$taxSpendSelf$", "$jobsForAll$", "$riskUnemployment$", "$riskPoverty$"),
                     no.space = TRUE,
                     omit.stat = c("all"),
                     covariate.labels = c("$LostEmployment$", "$UniversalCredit$", "$Furlough$", "$SelfEmployment$", "$LostEmployment \\cdot UC$"),
                     float = FALSE,
                     add.lines = list(c("Ind. FEs", "Yes", "Yes", "Yes", "Yes", "Yes"),
                                      c("Wave FEs", "Yes", "Yes", "Yes", "Yes", "Yes")))
sink()


## #################
## Descriptives

# Outcome presence/absence by wave

outcome_vars <- c("redistSelf", "taxSpendSelf", "ideoBHPSJjobForAll", "riskUnemployment", "riskPoverty", "govtHandouts", "reasonForUnemployment")

mat <- matrix(NA, length(outcome_vars), length(unique(bes_long$wave)))

for(i in 1:length(outcome_vars)) mat[i,] <- !sapply(unique(bes_long$wave),function(x) all(is.na(bes_long[[outcome_vars[i]]][bes_long$wave == x])))

mat <- data.frame(mat)

rownames(mat) <- outcome_vars
names(mat) <- unique(bes_long$wave)
mat$id <- outcome_vars
mat <- melt(mat, id.var = "id")

mat <- data.table(mat)
x <- mat[,sum(value),by = id]

mat$id[mat$id == "ideoBHPSJjobForAll"] <- "jobsForAll"
x$id[x$id == "ideoBHPSJjobForAll"] <- "jobsForAll"
setkey(x, V1)

mat$id <- factor(mat$id, levels = x$id)

mat$value <- ifelse(mat$value, "...in wave", "...not in wave")

dv_by_wave <- ggplot(mat, aes(x = variable, y = id, col = value, fill = value)) + geom_tile() + 
  scale_color_manual("Question...",values = c(scales::alpha("black",.7), scales::alpha("black",.2))) +
  scale_fill_manual("Question...",values = c(scales::alpha("black",.7), scales::alpha("black",.2))) + 
  theme_bw() + 
  theme(panel.grid = element_blank(), axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + 
  scale_x_discrete(labels = c(paste0("BES ",1:19),"PACER 1", "BES 20","PACER 2","BES 21")) + 
  xlab("") + ylab("")

png("../generated_images/DV_by_wave.png", 2300, 1000, res=  300)
print(dv_by_wave)  
dev.off()


# Raw outcome by wave

wave_dates <- as.Date(c("2014-02-20",
                        "2014-05-22",
                        "2014-09-19",
                        "2015-03-04",
                        "2015-03-31",
                        "2015-05-08",
                        "2016-04-14",
                        "2016-05-06",
                        "2016-06-24",
                        "2016-11-24",
                        "2017-04-24",
                        "2017-05-05",
                        "2017-06-09",
                        "2018-05-04",
                        "2019-03-11",
                        "2019-05-24",
                        "2019-11-01",
                        "2019-11-13",
                        "2019-12-13",
                        "2020-04-18",
                        "2020-06-03",
                        "2020-09-15",
                        "2021-05-08")) 

wave_dates <- data.frame(wave = c(1:19,19.5, 20, 20.5, 21), wave_dates)

bes_long <- merge(bes_long, wave_dates, by = "wave")
out <- merge(out, wave_dates, by = "wave")

plot_by_wave <- function(var = "taxSpendSelf", title = "Support for redistribution by wave"){
  
  title <- gsub("ideoBHPSJjobForAll", "jobsForAll", title)
  
  for_plot <- bes_long %>%
    mutate(outcome = as.numeric(get(var))) %>%
    filter(!is.na(cross_sectional_weight_) & !is.na(outcome)) %>%
    group_by(wave_dates) %>%
    mutate(n = sum(!is.na(outcome))) %>%
    filter(n != 1) %>%
    ungroup() %>%
    estimatr::lm_robust(outcome ~ as.factor(wave_dates) - 1, weights = cross_sectional_weight_, data = .) %>%
    DeclareDesign::tidy() %>%
    mutate(date = as.Date(gsub("as.factor\\(wave_dates\\)","",term))) %>% 
    filter(statistic != Inf) %>%
    mutate(estimate_mean = mean(estimate),
           estimate_sd = sd(estimate))
  
  for_plot %>%
    ggplot(aes(x = date, y = estimate, ymin = conf.low, ymax = conf.high)) + 
    geom_linerange() + 
    geom_line()  + 
    geom_point() + 
    theme_bw() + 
    theme(legend.position = "bottom",panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
    scale_x_continuous(breaks = as.Date(paste0(2014:2021,"-01-01")), labels = 2014:2021, limits = range(c(bes_long$wave_dates, as.Date("2021-01-01"))))  +
    geom_vline(xintercept = as.Date("2020-03-16"), lty = 2)+
    scale_linetype("", guide = "none") + 
    scale_size("", guide = "none") + 
    ggtitle(title) + xlab("Wave") + ylab("") 
  
}

plot_by_wave("redistSelf", "redistSelf") + ylim(c(6,7))
ggsave(paste0("../generated_images/fig1_redistSelf_by_wave.pdf"), width = 6, height = 3)
ggsave(paste0("../generated_images/fig1_redistSelf_by_wave.tiff"), width = 6, height = 3, dpi = 800)

plot_by_wave("taxSpendSelf", "taxSpendSelf") + ylim(c(6.75, 7.75))
ggsave(paste0("../generated_images/fig1_taxSpendSelf_by_wave.pdf"), width = 6, height = 3)
ggsave(paste0("../generated_images/fig1_taxSpendSelf_by_wave.tiff"), width = 6, height = 3, dpi = 800)

plot_by_wave("ideoBHPSJjobForAll", "ideoBHPSJjobForAll") + ylim(c(2.5, 3.5))
ggsave(paste0("../generated_images/ideoBHPSJjobForAll_by_wave.pdf"), width = 6, height = 3)

plot_by_wave("riskUnemployment", "riskUnemployment") + ylim(c(1.75, 2.75))
ggsave(paste0("../generated_images/fig4_riskUnemployment_by_wave.pdf"), width = 6, height = 3)
ggsave(paste0("../generated_images/fig4_riskUnemployment_by_wave.tiff"), width = 6, height = 3, dpi = 800)

plot_by_wave("riskPoverty", "riskPoverty") + ylim(c(2, 3))
ggsave(paste0("../generated_images/fig4_riskPoverty_by_wave.pdf"), width = 6, height = 3)
ggsave(paste0("../generated_images/fig4_riskPoverty_by_wave.tiff"), width = 6, height = 3, dpi = 800)

plot_by_wave("reasonForUnemployment", "reasonForUnemployment") + ylim(c(2.75, 3.75))
ggsave(paste0("../generated_images/fig6_reasonForUnemployment_by_wave.pdf"), width = 6, height = 3)
ggsave(paste0("../generated_images/fig6_reasonForUnemployment_by_wave.tiff"), width = 6, height = 3, dpi = 800)

plot_by_wave("govtHandouts", "govtHandouts") + ylim(c(2, 3))
ggsave(paste0("../generated_images/fig6_govtHandouts_by_wave.pdf"), width = 6, height = 3)
ggsave(paste0("../generated_images/fig6_govtHandouts_by_wave.tiff"), width = 6, height = 3, dpi = 800)


plot_by_wave_interaction <- function(var = "taxSpendSelf", interaction ="", title = "Support for redistribution by wave"){
  
  title <- gsub("ideoBHPSJjobForAll", "jobsForAll", title)
  
   for_plot <- out %>%
    mutate(outcome = as.numeric(get(var)),
           interaction = get(interaction)) %>%
    filter(!is.na(cross_sectional_weight) & !is.na(outcome)) %>%
    group_by(wave_dates) %>%
    mutate(n = sum(!is.na(outcome))) %>%
    filter(n != 1) %>%
    ungroup() %>%
    filter(!is.na(interaction)) %>%
    nest(data = -interaction) %>%
    mutate(mod = map(data, ~ estimatr::lm_robust(.$outcome ~ as.factor(.$wave_dates) - 1, weights = .$cross_sectional_weight)),
           tidied = map(mod, DeclareDesign::tidy)) %>%
    unnest(tidied) %>%
    mutate(date = as.Date(substring(term,24, 1000))) %>% 
    filter(statistic != Inf) 
  
  for_plot %>%
    ggplot(aes(x = date, y = estimate, ymin = conf.low, ymax = conf.high, col = interaction)) + 
    geom_linerange() + 
    geom_line()  + 
    geom_point() + 
    theme_bw() + 
    theme(legend.position = "bottom",panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
    scale_x_continuous(breaks = as.Date(paste0(2014:2021,"-01-01")), labels = 2014:2021, limits = range(c(out$wave_dates, as.Date("2021-01-01")))) + 
    geom_vline(xintercept = as.Date("2020-03-16"), lty = 2)+
    ggtitle(title) + xlab("Wave") + ylab("") +
    scale_color_manual("", values = c("orange", "darkgreen", "brown")[1:length(unique(for_plot$interaction))])
  
}


plot_by_wave_interaction("taxSpendSelf", "taxSpendSelf_prior_group", "taxSpendSelf by prior attitudes")
ggsave(paste0("../generated_images/taxSpendSelf_by_prior_by_wave.pdf"), width = 6, height = 4)

plot_by_wave_interaction("taxSpendSelf", "vote19", "taxSpendSelf by past vote")
ggsave(paste0("../generated_images/taxSpendSelf_by_past_vote_by_wave.pdf"), width = 6, height = 4)

plot_by_wave_interaction("redistSelf", "redistSelf_prior_group", "redistSelf by prior attitudes")
ggsave(paste0("../generated_images/redistSelf_by_prior_by_wave.pdf"), width = 6, height = 4)

plot_by_wave_interaction("redistSelf", "vote19", "redistSelf by past vote")
ggsave(paste0("../generated_images/redistSelf_by_past_vote_by_wave.pdf"), width = 6, height = 4)


## Descriptives of treatments
library(xtable)
out %>%
  filter(wave > 19) %>%
  group_by(wave) %>%
  rename(Wave = wave) %>%
  summarise(`Gov Support` = round(weighted.mean(income_guarantee_treatment, cross_sectional_weight) * 100, 2),
            `Furlough` = round(weighted.mean(income_guarantee_covid_job_treatment, cross_sectional_weight) * 100, 2),
            `Self-emp.` = round(weighted.mean(income_guarantee_covid_self_treatment, cross_sectional_weight) * 100, 2),
            `UC`  = round(weighted.mean(income_guarantee_other_treatment, cross_sectional_weight) * 100,2),
            `Lost Emp.` = round(weighted.mean(lost_employment, cross_sectional_weight) * 100,2),
            N = length(lost_employment)) %>%
  mutate(Wave = c("PACER 1 (April 2020)", "BES 20 (June 2020)", "PACER 2 (Sept 2020)")) %>%
  xtable(digits= 0,  caption = "Treatment distribution by wave (percentage of respondents per wave)", label = "tab:treat_dist_by_wave") %>%
  print.xtable(include.rownames = FALSE) %>%
  writeLines(con = "../generated_tex/treatment_variation.tex")


## Plot of treatment availability

outcome_vars <- c("income_guarantee_covid_binary", "income_guarantee_covid_binary", "income_guarantee_other_binary", "income_guarantee_binary", "workingStatus")

mat <- matrix(NA, length(outcome_vars), length(unique(bes_long$wave)))

for(i in 1:length(outcome_vars)) mat[i,] <- !sapply(unique(bes_long$wave),function(x) all(is.na(bes_long[[outcome_vars[i]]][bes_long$wave == x])))

mat <- data.frame(mat)

rownames(mat) <- c("CovidFurlough", "CovidSelfEmployment", "CovidUniversalCredit", "CovidGovSupport", "CovidLostEmployment")#outcome_vars
names(mat) <- unique(bes_long$wave)
mat$id <- c("CovidFurlough", "CovidSelfEmployment", "CovidUniversalCredit", "CovidGovSupport", "CovidLostEmployment")#outcome_vars
mat <- melt(mat, id.var = "id")

mat <- data.table(mat)
x <- mat[,sum(value),by = id]

mat$id[mat$id == "ideoBHPSJjobForAll"] <- "jobsForAll"
x$id[x$id == "ideoBHPSJjobForAll"] <- "jobsForAll"
setkey(x, V1)

mat$id <- factor(mat$id, levels = x$id)

mat$value <- ifelse(mat$value, "...in wave", "...not in wave")

iv_by_wave <- ggplot(mat, aes(x = variable, y = id, col = value, fill = value)) + geom_tile() + 
  scale_color_manual("Question...",values = c(alpha("black",.7), alpha("black",.2))) +
  scale_fill_manual("Question...",values = c(alpha("black",.7), alpha("black",.2))) + 
  theme_bw() + 
  theme(panel.grid = element_blank(), axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + 
  scale_x_discrete(labels = c(paste0("BES ",1:19),"PACER 1", "BES 20","PACER 2","BES 21")) + 
  xlab("") + ylab("")

png("../generated_images/IV_by_wave.png", 3000, 2500/2, res=  300)
print(iv_by_wave)  
dev.off()


### Parallel trends analysis

wave_by_wave1 <- out %>%
  filter(wave > 14 & !is.na(taxSpendSelf)) %>%
  mutate(taxSpendSelf = scale(as.numeric(taxSpendSelf))) %>%
  mutate(redistSelf = scale(as.numeric(redistSelf))) %>%
  group_by(id) %>%
  mutate(income_guarantee_treatment_all_waves = any(income_guarantee_treatment),
         lost_employment_all_waves = any(lost_employment)) %>%
  ungroup() %>%
  nest(data = -wave) %>%
  mutate(fit_1 = map(data, ~ broom::tidy(lm(scale(as.numeric(taxSpendSelf)) ~ -1 + lost_employment_all_waves, weights = cross_sectional_weight, data= .x), conf.int = TRUE)),
         fit_2 = map(data, ~ broom::tidy(lm(scale(as.numeric(taxSpendSelf)) ~ -1 + income_guarantee_treatment_all_waves, weights = cross_sectional_weight, data= .x), conf.int = TRUE)))

wave_by_wave2 <- out %>%
  filter(wave > 14) %>%
  mutate(taxSpendSelf = scale(as.numeric(taxSpendSelf))) %>%
  mutate(redistSelf = scale(as.numeric(redistSelf))) %>%
  group_by(id) %>%
  mutate(income_guarantee_treatment_all_waves = any(income_guarantee_treatment),
         lost_employment_all_waves = any(lost_employment)) %>%
  ungroup() %>%
  nest(data = -wave) %>%
  mutate(fit_3 = map(data, ~ broom::tidy(lm(scale(as.numeric(redistSelf)) ~ -1 + lost_employment_all_waves, weights = cross_sectional_weight, data= .x), conf.int = TRUE)),
         fit_4 = map(data, ~ broom::tidy(lm(scale(as.numeric(redistSelf)) ~ -1 + income_guarantee_treatment_all_waves, weights = cross_sectional_weight, data= .x), conf.int = TRUE)))
         
a <- wave_by_wave1 %>%
  select(wave, fit_1) %>%
  unnest(fit_1) %>%
  mutate(treat_var = "CovidLostEmployment",
         outcome = "taxSpendSelf") 

b <- wave_by_wave1 %>%
  select(wave, fit_2) %>%
  unnest(fit_2) %>%
  mutate(treat_var = "CovidGovSupport",
         outcome = "taxSpendSelf") 

c <- wave_by_wave2 %>%
  select(wave, fit_3) %>%
  unnest(fit_3) %>%
  mutate(treat_var = "CovidLostEmployment",
         outcome = "redistSelf") 

d <- wave_by_wave2 %>%
  select(wave, fit_4) %>%
  unnest(fit_4) %>%
  mutate(treat_var = "CovidGovSupport",
         outcome = "redistSelf") 

combined <- bind_rows(a,b,c,d) %>%
  mutate(term = ifelse(grepl("TRUE", term), TRUE, FALSE))
  
combined %>%
  ggplot(aes(x = wave, y = estimate, ymin = conf.low, ymax = conf.high, col = term)) + geom_pointrange() + geom_line() + 
  facet_grid(treat_var~outcome) + 
  #ylim(c(-.7,.7)) + 
  theme_bw() + 
  xlab("Wave") + 
  scale_x_continuous(breaks = c(1:19,19.5,20,20.5,21), labels = c(paste0("BES",1:19),"PACER W1", "BES20", "PACER W2", "BES21")) + 
  theme(panel.grid.minor = element_blank(), axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + 
  geom_vline(xintercept = 19.25, linetype = 2) + 
  scale_color_manual("",values = c("black", "red"), labels = c("Control", "Treatment")) + 
  ylab("Mean response") + 
  theme(legend.position = "bottom")

ggsave("../generated_images/parallel_trends_redistSelf_taxSpendSelf.pdf", width = 9, height = 6)

### Parallel trends analysis -- decomposed by prior attitude groups

library(weights, warn.conflicts = FALSE)

plot_parallel_trend_interaction <- function(outcome = "taxSpendSelf", interaction = "taxSpendSelf_prior_group"){

  out %>%
    mutate(Y = scale(as.numeric(get(outcome))),
           interaction = get(interaction)) %>%
    filter(!is.na(interaction)) %>%
    group_by(id) %>%
    mutate(income_guarantee_treatment_all_waves = any(income_guarantee_covid_self_treatment)) %>%
    ungroup() %>%
    group_by(wave, interaction) %>%
    nest(data = -c(wave, interaction)) %>%
    summarise(ttest = map(data, ~wtd.t.test(x = .$Y[.$income_guarantee_treatment_all_waves], 
                                            y = .$Y[!.$income_guarantee_treatment_all_waves],
                                            weight = .$cross_sectional_weight[.$income_guarantee_treatment_all_waves],
                                            weighty = .$cross_sectional_weight[!.$income_guarantee_treatment_all_waves],
                                            samedata = FALSE)),
              diff = map(ttest, ~.$additional[1]),
              se = map(ttest, ~.$additional[4])) %>%
    unnest(-c(ttest)) %>%
    mutate(hi = diff + 1.96 * se,
           lo = diff - 1.96 * se,
           wave = as.numeric(wave))  %>%
    unnest() %>%
    filter(!is.na(diff)) %>%
    ggplot(aes(x = wave, y = diff, col = interaction, ymin = lo, ymax = hi)) +
    geom_pointrange() + geom_line() + facet_wrap(~interaction) + 
    geom_hline(yintercept = 0, lty = 3) + 
    theme_bw() + 
    geom_vline(xintercept = 19, lty = 2) + 
    xlab("BES wave") + 
    ylab("Furloughed - not-furloughed") + 
    scale_color_manual("", values = c("orange", "darkgreen", "brown")) + 
    theme(legend.position = "bottom") + 
    ggtitle(outcome) + 
    xlim(1,20.5)
  
    
}

plot_parallel_trend_interaction("taxSpendSelf", "taxSpendSelf_prior_group")
ggsave("../generated_images/parallel_trends_taxSpendSelf_by_prior.pdf", width = 9, height = 6)

plot_parallel_trend_interaction("redistSelf", "redistSelf_prior_group")
ggsave("../generated_images/parallel_trends_redistSelf_by_prior.pdf", width = 9, height = 6)


## Change in unemployment risk

bes_long <- tibble(bes_long)
levels(bes_long$p_gross_personal) <- substring(levels(bes_long$p_gross_personal),2, 10000)
levels(bes_long$p_gross_household)[14:15] <- "£100,000 and over"

levels(bes_long$p_gross_personal) <- gsub(" per year","", levels(bes_long$p_gross_personal))
levels(bes_long$p_gross_household) <- gsub(" per year","", levels(bes_long$p_gross_household))

tmp_household <- bes_long %>% 
  transmute(riskUnemployment = as.numeric(scale(as.numeric(riskUnemployment))),
            p_gross_household = p_gross_household,
            wave = wave,
            weight = cross_sectional_weight_) %>%
  filter(wave %in% c(17,20)) %>%
  filter(!p_gross_household %in% c("Don't know", "Prefer not to answer", "988")) %>%
  filter(!is.na(p_gross_household)) %>%
  nest_by(p_gross_household) %>%
  mutate(mod = list(lm(riskUnemployment ~ as.factor(wave), data = data, weights = weight))) %>%
  mutate(ests = list(broom::tidy(mod)),
         est = as.numeric(ests[2,2]),
         se = as.numeric(ests[2,3]),
         hi = est + 1.96*se,
         lo = est - 1.96*se)

tmp_personal <- bes_long %>% 
  transmute(riskUnemployment = as.numeric(scale(as.numeric(riskUnemployment))),
            p_gross_household = p_gross_personal,
            wave = wave,
            weight = cross_sectional_weight_) %>%
  filter(wave %in% c(17,20)) %>%
  filter(!p_gross_household %in% c("Don't know", "Prefer not to answer", "988")) %>%
  filter(!is.na(p_gross_household)) %>%
  nest_by(p_gross_household) %>%
  mutate(mod = list(lm(riskUnemployment ~ as.factor(wave), data = data, weights = weight))) %>%
  mutate(ests = list(broom::tidy(mod)),
         est = as.numeric(ests[2,2]),
         se = as.numeric(ests[2,3]),
         hi = est + 1.96*se,
         lo = est - 1.96*se)

tmp_personal$type <- "Personal Income"
tmp_household$type <- "Household Income"

bind_rows(tmp_personal, tmp_household) %>% 
  ggplot(aes(y = p_gross_household, x = est, xmin = lo, xmax = hi)) + 
  geom_pointrange() + 
  theme_bw() + 
  geom_vline(xintercept = 0, lty = 2) + 
  xlab("Mean difference in riskUnemployment, wave 17 to wave 20") + 
  ylab("Yearly Gross Income") + 
  facet_wrap(~type)
ggsave("../generated_images/change_in_riskunemployment_by_income.pdf", width = 9, height = 6)

## Attrition

bes_long <- tibble(bes_long)

x <- bes_long %>%
  filter(id_in_wave == TRUE) %>%
  group_by(id) %>%
  summarise(id_in_crisis_wave = any(id_in_any_pacer == TRUE|id_in_bes_w20 == TRUE),
            id_in_pre_crisis_wave = min(wave[id_in_wave == TRUE]) <= 19,
            furlough_support = any(income_guarantee_covid_job_binary),
            self_employ_support = any(income_guarantee_covid_self_binary),
            other_support = any(income_guarantee_other_binary),
            p_sexuality = unique(p_sexuality),
            p_ethnicity = unique(p_ethnicity),
            p_religion = unique(p_religion),
            p_edlevel = unique(p_edlevel),
            p_parent = unique(p_parent),
            p_socgrade = unique(p_socgrade),
            p_marital = unique(p_marital),
            p_housing = unique(p_housing),
            p_gross_personal = unique(p_gross_personal)) %>%
  mutate(any_support = coalesce(furlough_support,self_employ_support, other_support),
         any_support = ifelse(!is.na(any_support), any_support, FALSE))

model_formula <- as.formula(any_support ~ p_gross_personal +  p_edlevel + p_socgrade  + p_ethnicity  + p_marital + p_housing)

pscore_model <- lm(any_support ~ p_gross_personal +  p_edlevel + p_socgrade , data=x, subset = x$id_in_pre_crisis_wave & x$id_in_crisis_wave)

x$pscore <- predict(pscore_model, newdata = x, type = "response")

#x %>%
#  group_by(any_support) %>%
#  summarise(mean(pscore, na.rm = TRUE))

# summary(lm(id_in_crisis_wave ~ pscore, x, subset = x$id_in_pre_crisis_wave))

## Work out the between-wave attrition for the most recent pre-crisis wave to the first post-crisis wave and compare that to the attirtion between pairs of other pre-crisis waves

wave_dates <- as.Date(c("2014-02-20",
                        "2014-05-22",
                        "2014-09-19",
                        "2015-03-04",
                        "2015-03-31",
                        "2015-05-08",
                        "2016-04-14",
                        "2016-05-06",
                        "2016-06-24",
                        "2016-11-24",
                        "2017-04-24",
                        "2017-05-05",
                        "2017-06-09",
                        "2018-05-04",
                        "2019-03-11",
                        "2019-05-24",
                        "2019-11-01",
                        "2019-11-13",
                        "2019-12-13",
                        #"2020-04-18",
                        "2020-06-03",
                        #"2020-09-15",
                        "2021-05-08"))

out_list <- list()
waves <- 1:19
for(w in waves){
  n_from_prior_wave <- length(bes_long$id[bes_long$wave == w & bes_long$id_in_wave == TRUE])
  n_from_prior_wave_in_this_wave <- sum(bes_long$id[bes_long$wave == w & bes_long$id_in_wave == TRUE] %in% bes_long$id[bes_long$wave == w+1 & bes_long$id_in_wave == TRUE])
  retention <- mean(bes_long$id[bes_long$wave == w & bes_long$id_in_wave == TRUE] %in% bes_long$id[bes_long$wave == w+1 & bes_long$id_in_wave == TRUE])  
  out_list[[w]] <- data.frame(date = wave_dates[w+1], 
                              retention = retention*100, 
                              wave = w+1, 
                              n_from_prior_wave, 
                              n_from_prior_wave_in_this_wave)
}


retention <- bind_rows(out_list) 

retention %>%
  ggplot(aes(x = date, y = retention)) + 
  geom_point() + geom_line() + 
  theme_bw() + 
  xlab("BES Wave") + 
  ylab("Retention from previous wave (%)") + 
  scale_x_date(breaks = retention$date, labels = 2:20) + 
  theme(panel.grid.major  = element_blank()) +
  geom_vline(xintercept = as.Date("2020-03-16"), lty = 2) +
  ylim(0,100)

ggsave("../generated_images/aggregate_retention.pdf", height = 4, width = 6)


## Distribution of government support by income group

x <- out %>% 
  filter(wave == 20 & id_in_wave == TRUE) %>%
  group_by(id) %>%
  summarise(received_support_furlough = any(income_guarantee_covid_treatment),
            received_support_self = any(income_guarantee_covid_self_treatment),
            received_support_UC = any(income_guarantee_other_treatment),
            p_gross_household = unique(p_gross_household),
            weight = unique(cross_sectional_weight)) %>%
  nest_by(p_gross_household) %>%
  mutate(mod_furlough = list(lm(received_support_furlough ~ 1, data = data, weights = weight)),
         ests_furlough = list(broom::tidy(mod_furlough)),
         est_furlough = as.numeric(ests_furlough[1,2]),
         se_furlough = as.numeric(ests_furlough[1,3]),
         hi_furlough = est_furlough + 1.96*se_furlough,
         lo_furlough = est_furlough - 1.96*se_furlough,
         
         mod_self = list(lm(received_support_self ~ 1, data = data, weights = weight)),
         ests_self = list(broom::tidy(mod_self)),
         est_self = as.numeric(ests_self[1,2]),
         se_self = as.numeric(ests_self[1,3]),
         hi_self = est_self + 1.96*se_self,
         lo_self = est_self - 1.96*se_self,
         
         mod_UC = list(lm(received_support_UC ~ 1, data = data, weights = weight)),
         ests_UC = list(broom::tidy(mod_UC)),
         est_UC = as.numeric(ests_UC[1,2]),
         se_UC = as.numeric(ests_UC[1,3]),
         hi_UC = est_UC + 1.96*se_UC,
         lo_UC = est_UC - 1.96*se_UC) %>%
  select(-c("data", starts_with("mod_"), starts_with("ests_")))

x %>%
  pivot_longer(-1) %>%
  mutate(type = case_when(grepl("est_", name) ~ "est",
                          grepl("se_", name) ~ "se",
                          grepl("hi_", name) ~ "hi",
                          grepl("lo_", name) ~ "lo"),
         var = case_when(grepl("_furlough", name) ~ "Furlough",
                         grepl("_UC", name) ~ "Universal Credit",
                         grepl("_self", name) ~ "Self-Employment")) %>%
  select(-name) %>%
  filter(!p_gross_household %in% c("Don't know", "Prefer not to answer", "988")) %>%
  filter(!is.na(p_gross_household)) %>%
  pivot_wider(id_cols = c(p_gross_household, var),
              names_from = type, values_from = value) %>%
  ggplot(aes(x = est, xmin = lo, xmax = hi, y = p_gross_household)) + 
  geom_pointrange() + 
  theme_bw() + 
  facet_wrap(~var) + 
  xlab("Fraction of respondents receiving support") + 
  ylab("")

ggsave("../generated_images/government_support_by_income.pdf", height = 4, width = 7)
